################################################################################
########################     Analysis and Plots     ############################
################################################################################

library(Matching)
library(tidyverse)

################################################################################

# Input is virginclip_pop.RData
# I omit the command to read in the dataset as it will need a directory.

rm(list=ls())
load("~/data/virginclip_pop.RData")

################################################################################

# 1. No trimming on forest cover in 2005

virginclip_pop <- virginclip_pop %>%
  filter(peat_depth_GW==0)

# outcome vars

virginclip_pop$did2015 <- virginclip_pop$forestha_2015 - virginclip_pop$forestha_2011
# virginclip_pop$did2006 <- virginclip_pop$forestha_2010 - virginclip_pop$forestha_2006

virginclip_pop$did2014 <- virginclip_pop$forestha_2014 - virginclip_pop$forestha_2011
# virginclip_pop$did2007 <- virginclip_pop$forestha_2010 - virginclip_pop$forestha_2007

virginclip_pop$did2013 <- virginclip_pop$forestha_2013 - virginclip_pop$forestha_2011
# virginclip_pop$did2008 <- virginclip_pop$forestha_2010 - virginclip_pop$forestha_2008

virginclip_pop$did2012 <- virginclip_pop$forestha_2012 - virginclip_pop$forestha_2011
# virginclip_pop$did2009 <- virginclip_pop$forestha_2010 - virginclip_pop$forestha_2009

virginclip_pop$didpre = virginclip_pop$forestha_2010 - virginclip_pop$forestha_2004

t1 = 2010
t11 = 2011
t0 = 2004

t2 = 2015
t3 = 2014
t4 = 2013
t5 = 2012


virginclip_pop$ddd2015 <- virginclip_pop$did2015 - ((t2-t11)/(t1-t0))*virginclip_pop$didpre
virginclip_pop$ddd2014 <- virginclip_pop$did2014 - ((t3-t11)/(t1-t0))*virginclip_pop$didpre
virginclip_pop$ddd2013 <- virginclip_pop$did2013 - ((t4-t11)/(t1-t0))*virginclip_pop$didpre
virginclip_pop$ddd2012 <- virginclip_pop$did2012 - ((t5-t11)/(t1-t0))*virginclip_pop$didpre

# ps model
ps_virginclip_pop <- glm(moratorium ~ forestha_2005 + forestha_2006 + forestha_2007 + 
                       forestha_2008 + forestha_2009 + forestha_2010 +  elev + slope + 
                       distroads + dist_cities + dist_palms_km + dist_timber_km + dist_logging_km +
                       agc_bcc_mean +  pop_2005 + pop_2010, family = binomial, data = virginclip_pop)
summary(ps_virginclip_pop)

virginclip_pop$pscore<- ps_virginclip_pop$fitted

# DDD models

PSM_DDD_virginclip_pop2015 <- Match( Y = virginclip_pop$ddd2015, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.001)

summary(PSM_DDD_virginclip_pop2015)

PSM_DDD_virginclip_pop2014 <- Match( Y = virginclip_pop$ddd2014, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.001)

summary(PSM_DDD_virginclip_pop2014)


PSM_DDD_virginclip_pop2013 <- Match( Y = virginclip_pop$ddd2013, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.001)

summary(PSM_DDD_virginclip_pop2013)

PSM_DDD_virginclip_pop2012 <- Match( Y = virginclip_pop$ddd2012, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.001)

summary(PSM_DDD_virginclip_pop2012)

rm(ps_virginclip_pop)

save.image("~/results/dryland-np-ddd-1215-0001.RData")


### peatland ###

rm(list=ls())
load("~/data/virginclip_pop.RData")

################################################################################

# 1. No trimming on forest cover in 2005

virginclip_pop <- virginclip_pop %>%
  filter(peat_depth_GW>0)

# outcome vars

virginclip_pop$did2015 <- virginclip_pop$forestha_2015 - virginclip_pop$forestha_2011
#virginclip_pop$did2006 <- virginclip_pop$forestha_2010 - virginclip_pop$forestha_2006

virginclip_pop$did2014 <- virginclip_pop$forestha_2014 - virginclip_pop$forestha_2011
#virginclip_pop$did2007 <- virginclip_pop$forestha_2010 - virginclip_pop$forestha_2007

virginclip_pop$did2013 <- virginclip_pop$forestha_2013 - virginclip_pop$forestha_2011
#virginclip_pop$did2008 <- virginclip_pop$forestha_2010 - virginclip_pop$forestha_2008

virginclip_pop$did2012 <- virginclip_pop$forestha_2012 - virginclip_pop$forestha_2011
#virginclip_pop$did2009 <- virginclip_pop$forestha_2010 - virginclip_pop$forestha_2009

virginclip_pop$didpre = virginclip_pop$forestha_2010 - virginclip_pop$forestha_2004

t1 = 2010
t11 = 2011
t0 = 2004

t2 = 2015
t3 = 2014
t4 = 2013
t5 = 2012


virginclip_pop$ddd2015 <- virginclip_pop$did2015 - ((t2-t11)/(t1-t0))*virginclip_pop$didpre
virginclip_pop$ddd2014 <- virginclip_pop$did2014 - ((t3-t11)/(t1-t0))*virginclip_pop$didpre
virginclip_pop$ddd2013 <- virginclip_pop$did2013 - ((t4-t11)/(t1-t0))*virginclip_pop$didpre
virginclip_pop$ddd2012 <- virginclip_pop$did2012 - ((t5-t11)/(t1-t0))*virginclip_pop$didpre

# ps model
ps_virginclip_pop <- glm(moratorium ~ forestha_2005 + forestha_2006 + forestha_2007 + 
                       forestha_2008 + forestha_2009 + forestha_2010 +  elev + slope + 
                       distroads + dist_cities + dist_palms_km + dist_timber_km + dist_logging_km +
                       agc_bcc_mean +  pop_2005 + pop_2010, family = binomial, data = virginclip_pop)
summary(ps_virginclip_pop)

virginclip_pop$pscore<- ps_virginclip_pop$fitted

# DDD models


PSM_DDD_virginclip_pop2015 <- Match( Y = virginclip_pop$ddd2015, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.001)

summary(PSM_DDD_virginclip_pop2015)

PSM_DDD_virginclip_pop2014 <- Match( Y = virginclip_pop$ddd2014, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.001)

summary(PSM_DDD_virginclip_pop2014)


PSM_DDD_virginclip_pop2013 <- Match( Y = virginclip_pop$ddd2013, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.001)

summary(PSM_DDD_virginclip_pop2013)

PSM_DDD_virginclip_pop2012 <- Match( Y = virginclip_pop$ddd2012, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.001)

summary(PSM_DDD_virginclip_pop2012)


# remove data and save results

rm(virginclip_pop, ps_virginclip_pop)
save.image("~/results/peatland-np-ddd-1215-0001.RData")